!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!    This file is part of ICTP RegCM.
!
!    Use of this source code is governed by an MIT-style license that can
!    be found in the LICENSE file or at
!
!         https://opensource.org/licenses/MIT.
!
!    ICTP RegCM is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!
!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

module mod_cmip6_cnrm

  use mod_intkinds
  use mod_realkinds
  use mod_message
  use mod_date
  use mod_stdio
  use mod_memutil
  use mod_kdinterp
  use mod_cmip6_helper
  use netcdf

  implicit none

  private

  character(len=*), parameter :: cnrm_version = 'v20181205'
  character(len=*), parameter :: cnrm_version1 = 'v20181206'
  character(len=*), parameter :: cnrm_version2 = 'v20191021'

  public :: read_3d_cnrm, read_2d_cnrm, read_fx_cnrm, read_sst_cnrm

  contains

    subroutine read_hcoord_cnrm(ncid,lon,lat)
      implicit none
      integer(ik4), intent(in) :: ncid
      real(rkx), pointer, contiguous, dimension(:), intent(inout) :: lon, lat
      integer(ik4) :: istatus, idimid, ivarid
      integer(ik4) :: nlon, nlat
      istatus = nf90_inq_dimid(ncid,'lon',idimid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lon dim')
      istatus = nf90_inquire_dimension(ncid,idimid,len=nlon)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire lon dim')
      istatus = nf90_inq_dimid(ncid,'lat',idimid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lat dim')
      istatus = nf90_inquire_dimension(ncid,idimid,len=nlat)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire lat dim')
      call getmem1d(lon,1,nlon,'cmip6_cnrm:lon')
      call getmem1d(lat,1,nlat,'cmip6_cnrm:lat')
      istatus = nf90_inq_varid(ncid,'lon',ivarid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lon var')
      istatus = nf90_get_var(ncid,ivarid,lon)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error read lon var')
      istatus = nf90_inq_varid(ncid,'lat',ivarid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lat var')
      istatus = nf90_get_var(ncid,ivarid,lat)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error read lat var')
    end subroutine read_hcoord_cnrm

    subroutine read_hcoord_sst_cnrm(ncid,lon,lat)
      implicit none
      integer(ik4), intent(in) :: ncid
      real(rkx), pointer, contiguous, dimension(:,:), intent(inout) :: lon, lat
      integer(ik4) :: istatus, idimid, ivarid
      integer(ik4) :: nlon, nlat
      istatus = nf90_inq_dimid(ncid,'x',idimid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lon x')
      istatus = nf90_inquire_dimension(ncid,idimid,len=nlon)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire x dim')
      istatus = nf90_inq_dimid(ncid,'y',idimid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find y dim')
      istatus = nf90_inquire_dimension(ncid,idimid,len=nlat)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire y dim')
      call getmem2d(lon,1,nlon,1,nlat,'cmip6_cnrm:lon')
      call getmem2d(lat,1,nlon,1,nlat,'cmip6_cnrm:lat')
      istatus = nf90_inq_varid(ncid,'lon',ivarid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lon var')
      istatus = nf90_get_var(ncid,ivarid,lon)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error read lon var')
      istatus = nf90_inq_varid(ncid,'lat',ivarid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lat var')
      istatus = nf90_get_var(ncid,ivarid,lat)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error read lat var')
    end subroutine read_hcoord_sst_cnrm

    subroutine read_vcoord_cnrm(ncid,ap,b)
      implicit none
      integer(ik4), intent(in) :: ncid
      real(rkx), pointer, contiguous, dimension(:), intent(inout) :: ap, b
      integer(ik4) :: istatus, idimid, ivarid
      integer(ik4) :: nlev
      istatus = nf90_inq_dimid(ncid,'lev',idimid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find lev dim')
      istatus = nf90_inquire_dimension(ncid,idimid,len=nlev)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire lev dim')
      call getmem1d(ap,1,nlev,'cmip6_cnrm:ap')
      call getmem1d(b,1,nlev,'cmip6_cnrm:b')
      istatus = nf90_inq_varid(ncid,'ap',ivarid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find a var')
      istatus = nf90_get_var(ncid,ivarid,ap)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error read a var')
      istatus = nf90_inq_varid(ncid,'b',ivarid)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error find b var')
      istatus = nf90_get_var(ncid,ivarid,b)
      call cmip6_error(istatus,__FILE__,__LINE__,'Error read b var')
    end subroutine read_vcoord_cnrm

    recursive subroutine read_3d_cnrm(idate,v,lonlyc)
      implicit none
      type(rcm_time_and_date), intent(in) :: idate
      type(cmip6_3d_var), pointer, intent(inout) :: v
      logical, optional, intent(in) :: lonlyc
      integer(ik4) :: istatus, idimid, it, irec
      integer(ik4) :: year, month, day, hour, y1, y2, m1, m2
      character(len=32) :: timecal, timeunit
      integer(ik4), dimension(4) :: istart, icount
      real(rk8), dimension(2) :: times
      type(rcm_time_interval) :: tdif
      character(len=16) :: ver

      if ( v%ncid == -1 ) then
        call split_idate(idate, year, month, day, hour)
        if ( v%vname == 'ua' .or. v%vname == 'va' ) then
          y1 = year
          m1 = ((month-1)/4) * 4 + 1
          if ( m1 == month .and. day == 1 .and. hour == 0 ) then
            m1 = m1 - 4
            if ( m1 < 1 ) then
              m1 = 12 + m1
              y1 = y1 - 1
            end if
          end if
          m2 = m1 + 4
          if ( m2 > 12 ) then
            m2 = 1
            y2 = y1 + 1
          else
            y2 = y1
          end if
        else
          y1 = year
          m1 = ((month-1)/6) * 6 + 1
          if ( m1 == month .and. day == 1 .and. hour == 0 ) then
            m1 = m1 - 6
            if ( m1 < 1 ) then
              m1 = 12 + m1
              y1 = y1 - 1
            end if
          end if
          m2 = m1 + 6
          if ( m2 > 12 ) then
            m2 = 1
            y2 = y1 + 1
          else
            y2 = y1
          end if
        end if
        if ( year < 2015 ) then
          ver = cnrm_version1
        else
          ver = cnrm_version2
        end if
        write(v%filename,'(a,i4,i0.2,a,i4,i0.2,a)') &
          trim(cmip6_path(year,'6hrLev',ver,v%vname)), &
          y1, m1, '010600-', y2, m2, '010000.nc'
#ifdef DEBUG
        write(stderr,*) 'Opening ',trim(v%filename)
#endif
        istatus = nf90_open(v%filename,nf90_nowrite,v%ncid)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error opening file '//trim(v%filename)//'.')
      end if
      if ( .not. associated(v%hcoord) ) then
        allocate(v%hcoord)
        call read_hcoord_cnrm(v%ncid,v%hcoord%lon1d,v%hcoord%lat1d)
      end if
      if ( .not. associated(v%vcoord) ) then
        allocate(v%vcoord)
        call read_vcoord_cnrm(v%ncid,v%vcoord%ak,v%vcoord%bk)
      end if
      if ( .not. associated(v%var) ) then
        v%ni = size(v%hcoord%lon1d)
        v%nj = size(v%hcoord%lat1d)
        v%nk = size(v%vcoord%ak)
        call getmem3d(v%var,1,v%ni,1,v%nj,1,v%nk,'cmip6:cnrm:'//trim(v%vname))
#ifdef DEBUG
        write(stderr,*) 'Input shape for ',trim(v%vname),' = ', &
          v%ni,'x',v%nj,'x',v%nk
#endif
        if ( present(lonlyc) ) then
          if ( lonlyc ) return
        end if
      end if
      if ( v%ivar == -1 ) then
        istatus = nf90_inq_varid(v%ncid,v%vname,v%ivar)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error searchong '//trim(v%vname)//' var in file '// &
          trim(v%filename)//'.')
      end if
      if ( v%nrec == -1 ) then
        istatus = nf90_inq_dimid(v%ncid,'time',idimid)
        call cmip6_error(istatus,__FILE__,__LINE__,'Error find time dim')
        istatus = nf90_inquire_dimension(v%ncid,idimid,len=v%nrec)
        call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire time dim')
        istatus = nf90_inq_varid(v%ncid,'time',it)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error searching time var in file '//trim(v%filename)//'.')
        istatus = nf90_get_att(v%ncid,it,"calendar",timecal)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time attribute calendar in file '//&
          trim(v%filename)//'.')
        istatus = nf90_get_att(v%ncid,it,"units",timeunit)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time attribute units in file '//trim(v%filename)//'.')
        istart(1) = 1
        icount(1) = 2
        istatus = nf90_get_var(v%ncid,it,times,istart(1:1),icount(1:1))
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time from file '//trim(v%filename)//'.')
        v%first_date = timeval2date(times(1),timeunit,timecal)
      end if

      tdif = idate - v%first_date
      irec = nint(tohours(tdif)/6.0) + 1

      if ( irec > v%nrec ) then
        istatus = nf90_close(v%ncid)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error close file '//trim(v%filename)//'.')
        v%ncid = -1
        v%ivar = -1
        v%nrec = -1
        irec = 1
#ifdef DEBUG
        write(stderr, *) 'Closed file, switching to the next in series...'
#endif
        call read_3d_cnrm(idate,v)
      end if
      istart(1) = 1
      istart(2) = 1
      istart(3) = 1
      istart(4) = irec
      icount(1) = size(v%var,1)
      icount(2) = size(v%var,2)
      icount(3) = size(v%var,3)
      icount(4) = 1
      istatus = nf90_get_var(v%ncid,v%ivar,v%var,istart,icount)
      call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error read variable '//v%vname//' from '//trim(v%filename)//'.')
    end subroutine read_3d_cnrm

    recursive subroutine read_2d_cnrm(idate,v,lonlyc)
      implicit none
      type(rcm_time_and_date), intent(in) :: idate
      type(cmip6_2d_var), pointer, intent(inout) :: v
      logical, optional, intent(in) :: lonlyc
      integer(ik4) :: istatus, idimid, it, irec
      integer(ik4) :: year, month, day, hour
      character(len=32) :: timecal, timeunit
      integer(ik4), dimension(3) :: istart, icount
      real(rk8), dimension(2) :: times
      type(rcm_time_interval) :: tdif

      if ( v%ncid == -1 ) then
        call split_idate(idate, year, month, day, hour)
        if ( year < 2000 ) then
          write(v%filename,'(a,a)') &
            trim(cmip6_path(year,'6hrLev',cnrm_version1,v%vname)), &
            '195001010600-200001010000.nc'
        else if ( year < 2015 ) then
          if ( year == 2000 .and. month == 1 .and. &
               day == 1 .and. hour == 0 ) then
            write(v%filename,'(a,a)') &
              trim(cmip6_path(year,'6hrLev',cnrm_version1,v%vname)), &
              '195001010600-200001010000.nc'
          else
            write(v%filename,'(a,a)') &
              trim(cmip6_path(year,'6hrLev',cnrm_version1,v%vname)), &
              '200001010600-201501010000.nc'
          end if
        else if ( year < 2065 ) then
          if ( year == 2015 .and. month == 1 .and. &
               day == 1 .and. hour == 0 ) then
            write(v%filename,'(a,a)') &
              trim(cmip6_path(year,'6hrLev',cnrm_version1,v%vname)), &
              '200001010600-201501010000.nc'
          else
            write(v%filename,'(a,a)') &
              trim(cmip6_path(year,'6hrLev',cnrm_version2,v%vname)), &
              '201501010600-206501010000.nc'
          end if
        else
          if ( year == 2065 .and. month == 1 .and. &
               day == 1 .and. hour == 0 ) then
            write(v%filename,'(a,a)') &
              trim(cmip6_path(year,'6hrLev',cnrm_version2,v%vname)), &
              '201501010600-206501010000.nc'
          else
            write(v%filename,'(a,a)') &
              trim(cmip6_path(year,'6hrLev',cnrm_version2,v%vname)), &
              '206501010600-210101010000.nc'
          end if
        end if
#ifdef DEBUG
        write(stderr,*) 'Opening ',trim(v%filename)
#endif
        istatus = nf90_open(v%filename,nf90_nowrite,v%ncid)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error opening file '//trim(v%filename)//'.')
      end if
      if ( .not. associated(v%hcoord) ) then
        allocate(v%hcoord)
        call read_hcoord_cnrm(v%ncid,v%hcoord%lon1d,v%hcoord%lat1d)
      end if
      if ( .not. associated(v%var) ) then
        v%ni = size(v%hcoord%lon1d)
        v%nj = size(v%hcoord%lat1d)
        call getmem2d(v%var,1,v%ni,1,v%nj,'cmip6:cnrm:'//trim(v%vname))
#ifdef DEBUG
        write(stderr,*) 'Input shape for ',trim(v%vname),' = ',v%ni,'x',v%nj
#endif
        if ( present(lonlyc) ) then
          if ( lonlyc ) return
        end if
      end if
      if ( v%ivar == -1 ) then
        istatus = nf90_inq_varid(v%ncid,v%vname,v%ivar)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error searchong '//trim(v%vname)//' var in file '// &
          trim(v%filename)//'.')
      end if
      if ( v%nrec == -1 ) then
        istatus = nf90_inq_dimid(v%ncid,'time',idimid)
        call cmip6_error(istatus,__FILE__,__LINE__,'Error find time dim')
        istatus = nf90_inquire_dimension(v%ncid,idimid,len=v%nrec)
        call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire time dim')
        istatus = nf90_inq_varid(v%ncid,'time',it)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error searching time var in file '//trim(v%filename)//'.')
        istatus = nf90_get_att(v%ncid,it,"calendar",timecal)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time attribute calendar in file '//&
          trim(v%filename)//'.')
        istatus = nf90_get_att(v%ncid,it,"units",timeunit)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time attribute units in file '//trim(v%filename)//'.')
        istart(1) = 1
        icount(1) = 2
        istatus = nf90_get_var(v%ncid,it,times,istart(1:1),icount(1:1))
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time from file '//trim(v%filename)//'.')
        v%first_date = timeval2date(times(1),timeunit,timecal)
      end if

      tdif = idate - v%first_date
      irec = nint(tohours(tdif)/6.0) + 1

      if ( irec > v%nrec ) then
        istatus = nf90_close(v%ncid)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error close file '//trim(v%filename)//'.')
        v%ncid = -1
        v%ivar = -1
        v%nrec = -1
        irec = 1
#ifdef DEBUG
        write(stderr, *) 'Closed file, switching to the next in series...'
#endif
        call read_2d_cnrm(idate,v)
      end if
      istart(1) = 1
      istart(2) = 1
      istart(3) = irec
      icount(1) = size(v%var,1)
      icount(2) = size(v%var,2)
      icount(3) = 1
      istatus = nf90_get_var(v%ncid,v%ivar,v%var,istart,icount)
      call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error read variable '//v%vname//' from '//trim(v%filename)//'.')
    end subroutine read_2d_cnrm

    recursive subroutine read_fx_cnrm(v)
      implicit none
      type(cmip6_2d_var), pointer, intent(inout) :: v
      integer(ik4) :: istatus

      v%filename = trim(cmip6_fxpath(cnrm_version,v%vname))
#ifdef DEBUG
      write(stderr,*) 'Opening ',trim(v%filename)
#endif
      istatus = nf90_open(v%filename,nf90_nowrite,v%ncid)
      call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error opening file '//trim(v%filename)//'.')
      allocate(v%hcoord)
      call read_hcoord_cnrm(v%ncid,v%hcoord%lon1d,v%hcoord%lat1d)
      v%ni = size(v%hcoord%lon1d)
      v%nj = size(v%hcoord%lat1d)
      call getmem2d(v%var,1,v%ni,1,v%nj,'cmip6:cnrm:'//trim(v%vname))
#ifdef DEBUG
      write(stderr,*) 'Input shape for ',trim(v%vname),' = ',v%ni,'x',v%nj
#endif
      istatus = nf90_inq_varid(v%ncid,v%vname,v%ivar)
      call cmip6_error(istatus,__FILE__,__LINE__, &
        'Error searchong '//trim(v%vname)//' var in file '// &
        trim(v%filename)//'.')
      istatus = nf90_get_var(v%ncid,v%ivar,v%var)
      call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error read variable '//v%vname//' from '//trim(v%filename)//'.')
      istatus = nf90_close(v%ncid)
      call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error close file '//trim(v%filename)//'.')
    end subroutine read_fx_cnrm

    recursive subroutine read_sst_cnrm(idate,v,lat,lon)
      implicit none
      type(rcm_time_and_date), intent(in) :: idate
      type(cmip6_2d_var), intent(inout) :: v
      real(rkx), pointer, contiguous, dimension(:,:), intent(in) :: lat, lon
      integer(ik4) :: istatus, idimid, it, irec
      integer(ik4) :: year, month, day, hour
      character(len=32) :: timecal, timeunit
      integer(ik4), dimension(3) :: istart, icount
      real(rk8), dimension(2) :: times
      type(rcm_time_interval) :: tdif

      if ( v%ncid == -1 ) then
        call split_idate(idate, year, month, day, hour)
        if ( year < 1950 ) then
          write(v%filename,'(a,a)') &
            trim(cmip6_path(year,'Oday',cnrm_version1,v%vname)), &
            '18500101-19491231.nc'
        else if ( year >= 1950 .and. year < 2015 ) then
          write(v%filename,'(a,a)') &
            trim(cmip6_path(year,'Oday',cnrm_version1,v%vname)), &
            '19500101-20141231.nc'
        else if ( year >= 2015 ) then
          write(v%filename,'(a,a)') &
            trim(cmip6_path(year,'Oday',cnrm_version2,v%vname)), &
            '20150101-21001231.nc'
        end if
#ifdef DEBUG
        write(stderr,*) 'Opening ',trim(v%filename)
#endif
        istatus = nf90_open(trim(v%filename),nf90_nowrite,v%ncid)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error opening file '//trim(v%filename)//'.')
      end if
      if ( .not. associated(v%hcoord) ) then
        allocate(v%hcoord)
        call read_hcoord_sst_cnrm(v%ncid,v%hcoord%lon2d,v%hcoord%lat2d)
        call h_interpolator_create(v%hint(1), &
                                   v%hcoord%lat2d,v%hcoord%lon2d,lat,lon)
        call getmem2d(v%var,1,size(v%hcoord%lon2d,1), &
                            1,size(v%hcoord%lat2d,2), &
                            'cmip6_cnrm:'//trim(v%vname))
      end if
      if ( v%ivar == -1 ) then
        istatus = nf90_inq_varid(v%ncid,trim(v%vname),v%ivar)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error searchong '//trim(v%vname)// &
          ' var in file '//trim(v%filename)//'.')
      end if
      if ( v%nrec == -1 ) then
        istatus = nf90_inq_dimid(v%ncid,'time',idimid)
        call cmip6_error(istatus,__FILE__,__LINE__,'Error find time dim')
        istatus = nf90_inquire_dimension(v%ncid,idimid,len=v%nrec)
        call cmip6_error(istatus,__FILE__,__LINE__,'Error inquire time dim')
        istatus = nf90_inq_varid(v%ncid,'time',it)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error searching time var in file '//trim(v%filename)//'.')
        istatus = nf90_get_att(v%ncid,it,"calendar",timecal)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time attribute calendar in file '//&
          trim(v%filename)//'.')
        istatus = nf90_get_att(v%ncid,it,"units",timeunit)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time attribute units in file '//trim(v%filename)//'.')
        istart(1) = 1
        icount(1) = 2
        istatus = nf90_get_var(v%ncid,it,times,istart(1:1),icount(1:1))
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error reading time from file '//trim(v%filename)//'.')
        v%first_date = timeval2date(times(1),timeunit,timecal)
      end if

      tdif = idate - v%first_date
      irec = int((tohours(tdif)+12)/24) + 1

      if ( irec > v%nrec ) then
        istatus = nf90_close(v%ncid)
        call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error close file '//trim(v%filename)//'.')
        v%ncid = -1
        v%ivar = -1
        v%nrec = -1
        irec = 1
        call read_sst_cnrm(idate,v,lat,lon)
      end if
      istart(1) = 1
      istart(2) = 1
      istart(3) = irec
      icount(1) = size(v%var,1)
      icount(2) = size(v%var,2)
      icount(3) = 1
      istatus = nf90_get_var(v%ncid,v%ivar,v%var,istart,icount)
      call cmip6_error(istatus,__FILE__,__LINE__, &
          'Error read variable '//trim(v%vname)// &
          ' from '//trim(v%filename)//'.')
      where ( v%var > 0.9E+20_rkx )
        v%var = -9999.0_rkx
      else where
        v%var = v%var + 273.15_rkx
      end where
    end subroutine read_sst_cnrm

end module mod_cmip6_cnrm

! vim: tabstop=8 expandtab shiftwidth=2 softtabstop=2
